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What determines the static force chains in stressed granular media? 
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The determination of the normal and transverse (frictional) inter-particle forces within a granular 
medium is a long standing, daunting, and yet unresolved problem. We present a new formalism 
which employs the knowledge of the external forces and the orientations of contacts between particles 
(of any given sizes), to compute all the inter-particle forces. Having solved this problem we exemplify 
the efficacy of the formalism showing that the force chains in such systems are determined by an 
expansion in the eigenfunctions of a newly defined operator. 


In a highW influential paper from 2005 Majmudar and 
Behringer [l[ wrote: “Inter-particle forces in granular 
media form an inhomogeneous distribution of filamen¬ 
tary force chains. Understanding such forces and their 
spatial correlations, specifically in response to forces at 
the system boundaries, represents a fundamental goal of 
granular mechanics. The problem is of relevance to civil 
engineering, geophysics and physics, being important for 
the understanding of jamming, shear-induced yielding 
and mechanical response.” A visual example of such force 
chains in a system of plastic disks is provided in Fig. [T] 
In this Letter we present a solution of this goal. 

To be precise, the problem that we solve is the fol¬ 
lowing: consider a granular medium with known sizes of 
the granules, for example the 2-dimensional systems ana¬ 
lyzed in Ref. [H and shown in Fig. [1] of disks of known 
diameters {cri}A ^. Given the external forces, denoted be¬ 
low as and the external torques F™* exerted on the 
granules, and given the angular orientations of the vec¬ 
tors connecting the center of masses of contacting gran¬ 
ules (but not the distance between them!), determine all 
the inter-particle normal and tangential forces and 
/b. The method presented below applies to granular 



FIG. 1: Force chains in a binary system of plastic disks of 
two diameters, stressed uniaxially at the boundaries. The 
force chains are made visual by the optical birefringence of the 
stressed disks. The image is courtesy of V. Sathish Akella and 
Mahesh Bandi, Okinawa Institute of Science and Technology. 


systems in mechanical equilibrium; the issue of instabil¬ 
ities and abrupt changes in the force chains will be dis¬ 
cussed elsewhere. For the sake of clarity and simplicity 
we will present here the two-dimensional case; the savvy 
reader will recognize that the formalism and the solution 
presented will go smoothly also for the three-dimensional 
case (as long as the system is in mechanical equilibrium). 
The full formalism will be presented in a longer publica¬ 
tion in due course. 

The obvious conditions for mechanical equilibrium are 
that the forces and the torques on each particle have 
to sum up to zero 0. The condition of force balance 
is usefully presented in matrix form using the following 
notation. Denote the (signed) amplitudes of the inter¬ 
particle forces fij as a vector |/), where the amplitudes 
/”■ appear first and then the amplitudes /b. The vector 
of X and y components and is denoted as 
where all the x components are presented in first 

and then all the y components. The vector |/) has 2c 
entries where c is the number of contacts between par¬ 
ticles. The vector has 2N entries where N is the 

number of particles, with zero entries for all the particles 
on which there is no external force. We can then write 
the force balance condition as 

M\f) = , (1) 

where M is a 2N x 2c matrix. The entries in the matrix 
M contain the directional information, see Supplemen¬ 
tal Material at [URL will be inserted by publisher] for 
an example of an M matrix. Denote the unit vector in 
the direction of the vector distance between the centers 
of mass of particles i and j by n^, and the tangential 
vector by tij orthogonal to n^. Then the entries of M 
display the projections fiij^x and fiij^y or and Uj^y as 
appropriate. We thus guarantee that Eq. CD is equivalent 
to the mechanical equilibrium condition 

j2f,y + Fr^ = o. ( 2 ) 

j 

As is well known, the friction-less granular system in the 
thermodynamic limit is jammed exactly at the isostatic 
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condition 2N = c [^. In the friction-less case M is a 
2N X c matrix and as long as c = 2N one can solve the 
problem by multiplying Eq. CD by the transpose 
getting 

M^M\f) = . (3) 

In this case the matrix M'^M has generically exactly 
three Goldstone modes (two for translation and one for 
rotation) Q, and since the external force vector is or¬ 
thogonal to the Goldstone modes (otherwise the external 
forces will translate or rotate the system), Eq. ([3]) can be 
inverted with impunity by multiplying by In 

fact even when c < 2N but the system is small enough 
to be jammed, this method can be used since there are 
enough constraints to solve for the forces. This last com¬ 
ment is important for our developments below. 

The problem becomes under-determined above iso- 
staticity in the frictionless case, when force chains begin 
to build up that span from one boundary to the other. 
With friction we anyway have twice as many unknowns 
and we need to add the constrains of torque balance. The 
condition of torque balance is Vi x /h -\- T®’'* = 0 on 
every particle, where is the external torque exerted 
on the ith disk Q. Eor disks, is in the normal direc¬ 
tion, and therefore the torque balance becomes a condi¬ 
tion that the sum of tangential forces has to balance the 
external tangential force. This condition can be added 
to Eq. ([1]) using a new matrix B in the form 



The order of the extended matrix B is 3iV x 2c, see Sup¬ 
plemental Material at [URL will be inserted by publisher] 
for an example of T. Above jamming when the number 
of contacts increases 2c 3> 3iV. The matrix B is not 
square, and the matrix B which is of size 2c x 2c, has 
at least 2c — 3N zero modes @ . Accordingly it cannot 
be inverted and one can conclude that the conditions 
of mechanical equilibrium are not sufficient to de¬ 
termine all the forces. 

Obviously what is missing are additional constraints 
to remove the host of zero modes. These additional con¬ 
straints are geometrical constraints in which can be 
read from those disks which describe connected polygons. 
In other words, since we know the orientation of each 
contact in our system, we can determine which granules 
are stressed in a triangular arrangement, and which in 
a square or pentagonal etc., see Eig. |21 Each such ar¬ 
rangement is a constraint on the radius vectors adjoining 
the centers of mass. For example if particles i, j and k 
are in a triangular arrangement then -f Vjk -f r^i = 0, 
with the analogous constraint on squares, pentagons etc. 



FIG. 2: A generic situation in a stressed arrangement of 
N = 242 binary disks of diameters ai = 1, CT 2 = 1-5. 
Here c = 432. Shown are the P = 205 polygons: 80 tri¬ 
angles in blue, 77 squares in red, 36 pentagons in green, 7 
hexagons in yellow, 4 heptagons in cyan and 1 octagon in 
black. The Euler characteristic of the system is as expected 
A — c-|-(P-l-l) = 2-|-P since we have 14 rattlers. 

These constraints can be written in a matrix form by de¬ 
noting the amplitudes of inter-particle vector distances as 
jr) where we again arrange the x components first and 
the y components second: 

Q\r) = 0 , (5) 

where the matrix Q again has entries or as ap¬ 
propriate to represent the vectorial geometric constraints, 
see Supplemental Material at [URL will be inserted by 
publisher] for an example of Q. Denoting the total num¬ 
ber of polygons by P the dimension of the matrix Q is 
2P X c. Of course jr) has c entries while [/) had 2c en¬ 
tries. Note that in generic situations there can be also 
disks which are not stressed at all. These are referred to 
as “rattlers”. For example in the configuration shown in 
Fig. [2] there exist 14 rattlers. 

At this point we specialize the treatment to Hookean 
normal forces with a given force constant k Q. Non 
Hookean forces result in a nonlinear theory that can still 
be solved but much less elegantly. For the present case 

= K[{ai + crj)h,jj2 - . ( 6 ) 

Denoting the amplitudes of the vectors (tTi -|- aj)nijj2 
as the vector jcr) (again with first the x and then the y 
components), we can rewrite Eq. (I5|) in the form 

Q\r) = Q\a-r/K)=0 , (7) 

Q\r) = . ( 8 ) 

Having this result at hand we can formulate the final 
problem to be solved. Arrange now a new matrix, say 
G, operating on a vector [/), with a RHS being a vector, 
say Jt), made of a stacking of jF®’'*), jT®’'*) and Q\kct), 
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as before with x and then y components: 


G\f) = 



M 



( /r* 

xext 

pext 




\Q\Ka)/ 


Using these objects our problem is now 


( 9 ) 


G\f) = \t) . (10) 


The dimension of the matrix G is {3N + 2P) x 2c and the 
matrix G^G has the dimension 2c x 2c. We can use now 
the Euler characteristic @ to show that the situation has 
been returned here to the analog of the invertible matrix 
when c < 27V: the Euler characteristic in two 
dimensions requires that 


7V-c+(P + l) = 2 + i?, (11) 


where R is the number of “rattlers” i.e. disks on which 
there is no force. Accordingly we find that 

2c = 27V + 2P-2-2i?<37V + 2P . (12) 

Consequently, the matrix G^G has no zero eigenmodes. 
Thus the final solution for the forces can be obtained as 

|/) = (G^G)-iG^|f) = ^^Hlg!^|vi/J . (13) 

i 

where is the set of eigenfunctions of G^G associ¬ 
ated with eigenvalues . We compared the inter-particle 
forces obtained from direct numerical simulations (see be¬ 
low for details) to those computed from Eq. (fOl) . Both 
normal and tangential forces are of course identical to 
machine accuracy. We reiterate that we did not need to 
know the distances between particles. This is important 
in applying the formalism to experiments since it is very 
difficult to measure with precision the degree of com¬ 
pression of hard particles like, say, metal balls or sand 
particles. Note also the remarkable fact that we never 
had to provide the frictional (tangential) force law in the 
formalism to obtain the correct forces! 

At this point we can discuss the force chains. By defini¬ 
tion these are the large forces in the system that provide 
the tenuous network that keeps the system rigid. Observ¬ 
ing Eq. (IT^ we should focus on the eigenfunction vki of 
G^G that have the smallest eigenvalues and the largest 
overlaps with G'^\t). These can be found and arranged in 
order of the magnitude of i\G'^\t)/\ independently of 
the calculation of |/). In Fig. [3] we show the contribution 
to the total energy (/!/)/«:, learning that about 20% of 
the leading eigenfunction are responsible for 90% of the 
energy. We can therefore hope that the force chains will 
be determined by the same relatively small number of 



FIG. 3: The cumulative percentage contribution to the energy 
of the eigenfunctions of G^G, ordered according to the 
magnitude of (^i|G^|t}/Ai. The convergence is relatively fast 
with the first 168 leading eigenfunction (out of 864 modes) 
contributing 90% of the total energy. 


eigenfunctions. This is not guaranteed; due to contribu¬ 
tions to the forces that oscillate in sign the convergence 
can be much slower than in the case of the energy where 
the sum is of positive contributions. In Fig.|T]we show in 
upper left panel the force chains in the configuration of 
Fig. [2] In the other panels we show the prediction of the 
force chains using 100, 200 and 300 of the (energy) lead¬ 
ing modes. We learn that with 100 out of the 864 modes 
the main force chains begin to be visible. With 200 out 
of the 864 modes the full structure of the force chains 
is already apparent, although with 300 it is represented 
even better. 


Forces as computed from simulation Forces computed from 100 leading modes 



-10 -5 0 5 10 -10 -5 0 5 10 

Forces computed from 200 leading modes Forces computed from 300 leading modes 




-10 -5 0 5 10 


-10 -5 0 5 10 



FIG. 4: Upper panel: the force chains in the configuration 
shown in Fig. [2] Lower panel: the force chains as predicted 
from 100, 200 and 300 (energy) leading modes. While the 
main contributions to the force chains are visible already with 
100 out of the 864 modes, the full structure is apparent only 
with 200. Using 300 modes we see a reasonably faithful re¬ 
production of the force chains. 


Since the number of geometric constraints is very large, 
one can ask whether all the geometric constraints are 
necessary, as Eq. (ITT]) shows that 2c <C 37V -|- 2P. The 
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answer is no, we could leave out constraints as long as we 
have enough conditions to determine the solution. There 
is the obvious question then why do we have a unique 
solution when the number of equations is larger than the 
number of unknowns. The answer to this question lies in 
the properties of the vector |t) and the matrix GG^ which 
does have many zero modes. A condition for the existence 
of a solution is that \t) is orthogonal to all the zero modes 
of GG^, as can be easily checked. We have ascertained 
in our simulations that this condition is always met. 

In the near future we will present an extension of this 
formalism to three dimensions and the use of the for¬ 
malism to study the instabilities of the force networks to 
changes in the external forces. As a final comment we 
should note that in fact only one external force is nec¬ 
essary to determine all the inter-disk forces. This single 
external force is necessary to remove the re-scaling free¬ 
dom that this problem has by definition. 

Simulations: For the numerical experiments in 2- 
dimensions we use disks of two diameters, a ‘small’ one 
with diameter ag = 1.0 and a ‘large’ one with diameter 
af: = 1.5. Such N disks are put between virtual walls at 
X = ±a and y = ±&. These walls exert external forces on 
the disks. The external forces are taken as Hookean for 
simplicity. For disks near the wall at x = ±a we write 

= -iri,x-a) 

= -{ri,x + a) iiri^x<-a 

= 0 otherwise . (14) 

Here denotes the x component of the position vector 
ri of the center of mass of the ith disk, and we have a 
similar equation for the y components with a replaced 
by b. When two disks, say disk i and disk j are pressed 
against each other we define their amount of compression 
as 5rij\ 

5nj = (Tij - r^j , (Ty = {ai -I- aj)/2 , (15) 

where is the actual distance between the centers of 
mass of the disks i and j. 

In our simulations the normal force between the disks 
acts along the radius vector connecting the centers of 
mass. We employ a Hookean force = nSrijfiij. 

To define the tangential force between the disks we 
consider (an imaginary) tangential spring at every con¬ 
tact which is put at rest whenever a contact between the 
two disks is formed. During the simulation we implement 
memory such that for each contact we store the signed 
distance Stij to the initial rest state. For small deviations 
we require a linear relationship between the displacement 
and the acting tangential force. This relationship breaks 
when the magnitude of the tangential force reaches 
where due to Coulomb’s law the tangential loading can no 
longer be stored and is thus dissipated. When this limit 


is reached the bond breaks and after a slipping event the 
bond is restored with a the tangential spring being loaded 
to its full capacity (equal to the Coulomb limit). 

The numerical results reported above were obtained 
by starting with N = 242 particles on a rectangular grid 
(ratio 1:2) with small random deviations in space and 
no contacts. We implement a large box that contains all 
the particles. The box acts on the system by exerting a 
restoring harmonic normal force as described in Eq. m- 
The experiment is an iterative process in which we first 
shrink the containing box infinitesimally (conserving the 
ratio). The second step is to annul all the forces and 
torques, to bring the system back to a state of mechan¬ 
ical equilibrium. We therefore annul the forces using a 
conjugate gradient minimizer acting to minimize the re¬ 
sulting forces and torque on all particles. We iterate these 
two steps until the system is compressed to the desired 
state. 

This work had been supported in part by an “ideas” 
grant STANPAS of the ERC. We thank Deepak Dhar for 
some very useful discussions. We are grateful to Edan 
Lerner for reading an early version of the manuscript 
with very useful remarks. 
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Supplementary Material: Examples of the M,T and 
Q matrices 

We consider a two-dimensional configuration of N par¬ 
ticles with C contacts and P polygons. For convenience 
of notation, only single digit particle indices are used in 
this example, so that the notation hi 3 ^ means the Carte¬ 
sian X component of the unit vector from the center of 
particle 1 to that of particle 3. 

The convention for ordering of the contacts is demon¬ 


strated in Eq. [1^] (and see also Fig. [5]): 


c= 1 2 ... 7 8 ... \ . . 

Iwith2 IwithS ... 2with3 2with4 ... J 



FIG. 5: A subset of the particle configuration of the system. 
The particles’ diameters in the figure are scaled-up to visu¬ 
ally stress their finite overlap. Some particle indices used in 
the M and T matrices are shown. Arrows represent the nor¬ 
mal vectors used to construct the M and T matrices (before 
normalization). Different arrow colors are for visualization 
purposes only. 


The M matrix is used to describe the force balance condi¬ 
tion (Eq. 1 in the main text) and has dimension 2N x 2C 
in the most general case when contact forces have both 
normal and tangential components. Each row is associ¬ 
ated with a given particle i and each column describes 
one contact and has non-zero entries corresponding only 
to the pair of particles i and j forming that contact. Its 
first N rows store the x components and the next N 
rows store the y components of unit normal vectors fzy 
and unit tangential vectors (counter-clockwise orthog¬ 
onal to fiij). The first C columns of M correspond to the 
normal directions and the next C columns correspond to 
the tangential directions (which can also of course be ex¬ 
pressed using the normal directions via a simple rotation 
transformation). An example of some of the terms of the 
M matrix for the configuration of Fig. [5] is given in Eq. 

[HI 


/ -hi2^ 

—ni3 

0 

0 

.. il2" 

- 

0 

0 


0 

A X 
— ^23 

A X 
— ^24 

.. -ii2^ 

0 ... 

CO 

i24“ 

0 

/A ^ 
ni3 

/A X 
^23 

0 

.. 0 

-tl3" - 

-i23" 

0 

0 

0 

0 

A X 

7224 

.. 0 

0 ... 

0 

-t24" 

—ni2^ 

—ni3^ ... 

0 

0 

ti2^ 

tl3^ ... 

0 

0 

hi2^ 

0 

—n23^ 

-7224^ ... 

.. —ii2'^ 

0 ... 

i23^ 

i24^ 

0 

ni3^ 

'^23^ 

0 

.. 0 

—fl3^ ... 

— i23^ 

0 

0 

0 

0 

7124^ 

.. 0 

0 ... 

0 

—^24^ 




The T matrix is used to describe the torque balance con- given in Eq. [THl 
dition (see Eq. 9 in the main text) and is of dimensions 
N X C. Again, the row indices correspond to particles 


and the column indices refer to contacts. The non-zero 


/ Ri 

Ri 

. 0 

0 

entries in each column correspond to the radii of particles 


i?2 

0 . 

. R2 

i?2 

i and j forming that contact. An example of some of the 

T = 

0 

i?3 . 

. R3 

0 

terms of the T matrix for the configuration of Fig. [5] is 


0 

0 . 

. 0 

i?4 


V : : : : / 


(17) 


( 18 ) 
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When the external torque is zero, as in our loading proto¬ 
col using compression, the radii are eliminated from the 
torque balance equation and the T matrix can be further 


simplified to the form of 

Eq. 

M 




(1 

1 . 

. 0 

0 . 



1 

0 . 

. 1 

1 . 


T = 

0 

1 . 

. 1 

0 . 



0 

0 . 

. 0 

1 . 







/ 


The Q matrix (cf. Eq. 7 in the main text) is used to de¬ 
scribe the presence of closed polygons formed by particles 
in contact and and is of dimensions 2P x C. Here row 
indices correspond to polygons and column indices refer 
to the contacts. Non-zero entries in each row describe 
the unit normal directions joining two particles in con¬ 
tact which are members of a given polygon. The first P 
rows store the x components and the next P rows store 
the y components of unit vectors n^. An example for 
some of the terms of the Q matrix is given in Eq. [SD] 
(and see Fig. (5]): 


/ 

0 

.A ^ 
—ni3 

fils'" 

0 

A X 
—ni5 

A a: 

• n23 

. 0 . 

. 0 

A X 
■ «35 


hi2^ 

0 

V ^ 

-fiisy 

nis^ 

0 

• fl23^ 

. 0 . 

. 0 

■ ^35^ 

/ 


( 20 ) 



FIG. 6: Same configuration as in Fig. [5l only now the poly¬ 
gons used in the Q matrix are demonstrated. 














